function [ minGraphV ] = ExtendMinGraph(vSignature,...
                                        minGraphV,fullGraphV,fullGraphVDistanceMtr,...
                                        userFtPt,vertSurfArea)
% Extend the minimal graph until the matching to the full graph becomes inambiguous
% The method uses adajcency matrices
% vSignature:       signature of the vertices
% distanceMtr:      matrix of geodesic distances
% minGraphVert:     index of the vertices used for the minimal graph
% minGraphAdjMtr:   adjacency matrix of the minimal graph
% fullGraphVert:    index of the vertices used for the full graph
% fullGraphAdjMtrF: adjacency matrix of the full graph



initMinGraphVert=minGraphV;

minGraphVDistanceMtr = CompMinGraphVDistanceMtr(minGraphV,fullGraphV,fullGraphVDistanceMtr);

% Comute all possible matching
allMatch = ComputeMatching(minGraphV,vSignature(minGraphV,:),minGraphVDistanceMtr(:,minGraphV),...
                           fullGraphV,vSignature(fullGraphV,:),fullGraphVDistanceMtr(:,fullGraphV),intmax);
% Compute the matching cost for all these possible matching
cost = MatchingCost(allMatch,minGraphV,vSignature(minGraphV,:),minGraphVDistanceMtr(:,minGraphV),...
                             fullGraphV,vSignature(fullGraphV,:),fullGraphVDistanceMtr(:,fullGraphV));

% Sort them by decreasing order of quality
[~,IX]=sort(cost(:,2),'descend');
allMatch = allMatch(IX,:);
cost=cost(IX,:);

% Remove the matching of lowest cost -> exact matching
allMatch=allMatch(1:numel(allMatch)-1);

% for i=1:size(allMatch,1)
%     minGraphVtoFullGraphV=cell2mat(allMatch(i,1));
%     str=strcat('Found one matching: ',num2str(minGraphVtoFullGraphV(1)));
%     for j=2:numel(minGraphVtoFullGraphV)
%         str=strcat(str,',',num2str(minGraphVtoFullGraphV(j)));
%     end
%     disp(str);
% end

markedFullV=false(numel(fullGraphV),1);
for fullV=1:numel(fullGraphV)
    if find(minGraphV==fullGraphV(fullV))
        markedFullV(fullV)=true;
    end
end

if isempty(allMatch)
    prevMinGraphV=minGraphV;
else
    prevMinGraphV=minGraphV;
    prevAllMatch=allMatch;

    while size(prevAllMatch,1)>0
        [newMinGraphV,markedFullV] = AddVToMinGraph(prevMinGraphV,fullGraphV,fullGraphVDistanceMtr(:,fullGraphV),markedFullV);
        if numel(newMinGraphV)>numel(prevMinGraphV)
            % Compute the number of matching
            newMinGraphVDistanceMtr = CompMinGraphVDistanceMtr(newMinGraphV,fullGraphV,fullGraphVDistanceMtr);
            allMatch = CompleteMatching(newMinGraphV,vSignature(newMinGraphV,:),newMinGraphVDistanceMtr(:,newMinGraphV),...
                                        fullGraphV,vSignature(fullGraphV,:),fullGraphVDistanceMtr(:,fullGraphV),...
                                        minGraphV,prevAllMatch);
                       
            % If the quality is imptoved, keep the newly added feature point
            if size(allMatch,1)<size(prevAllMatch,1)
                prevMinGraphV=newMinGraphV;
                prevAllMatch=allMatch;
            end
        else
            break;
        end;
    end
end
minGraphV=prevMinGraphV;
minGraphVDistanceMtr = CompMinGraphVDistanceMtr(minGraphV,fullGraphV,fullGraphVDistanceMtr);

%%%%%%%%%%%%%%%%%%%%%
% Save all the data %
%%%%%%%%%%%%%%%%%%%%%
srcMinGraphV=minGraphV;                                     % Min graph on the source mesh
srcMinGraphVSignature=vSignature(minGraphV,:);              % Vertex signature
srcMinGraphFullDistanceMtr=minGraphVDistanceMtr;            % Geodesic distance of the extended min graph to all other vertices
srcUserFtPtID=userFtPt;
% Vasyl {
%srcInitMinGraphV={initMinGraphVert};
srcInitMinGraphV={initMinGraphVert};
% avoid using cell data structure to maintain compliance with Octave
srcInitMinGraphV = initMinGraphVert; 
% } Vasyl
srcVertSurfArea=vertSurfArea;


% Binary saving
save('dataForTrg','srcMinGraphV','srcMinGraphVSignature',...
    'srcMinGraphFullDistanceMtr','srcUserFtPtID','srcInitMinGraphV','srcVertSurfArea', 'userFtPt');

end


% Compute all the possible matching of minGraph to fullGraph
% Return these matching in allMatch
function [ allMatch ]= ComputeMatching(minGraphV,minGraphVSignature,minGraphVDistanceMtr,...
                                       fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr,nbMaxMatching)

global matchVThreshold;

allMatch={};

% The first vertex to match
minV=1;

% Compute all possible matching of the first vertex
for fullV=1:numel(fullGraphV)
    % Check if the signature matches
    if CompVertSimilarity(minGraphVSignature(minV),fullGraphVSignature(fullV)) < matchVThreshold
        continue;
    end
    % assignement vector that shows which vertex of min graph is assigned to which vertex of full graph
    minGraphVtoFullGraphV=zeros(numel(minGraphV),1);
    minGraphVtoFullGraphV(minV)=fullV;
    fullGraphVtoMinGraphV=zeros(numel(fullGraphV),1);
    fullGraphVtoMinGraphV(fullV)=minV;
    [allMatch]=RecursComputeMatching(minGraphVtoFullGraphV,fullGraphVtoMinGraphV,...
                                     minGraphV,minGraphVSignature,minGraphVDistanceMtr,...
                                     fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr,...
                                     allMatch,nbMaxMatching);
    if size(allMatch,2) >= nbMaxMatching
        break;
    end
end
end

% Continue the computation of the matching 
% Compute all the possible matching of minGraph to fullGraph
% Return these matching in allMatch
function [ allMatch ]= CompleteMatching(minGraphV,minGraphVSignature,minGraphVDistanceMtr,...
                                        fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr,...
                                        prevMinGraphV,prevAllMatch)
allMatch={};

for i=1:size(prevAllMatch,1)
    
    % Put the already matching vertices
    minGraphVtoFullGraphV=zeros(numel(minGraphV),1);
    fullGraphVtoMinGraphV=zeros(numel(fullGraphV),1);
    prevFullGraphV=cell2mat(prevAllMatch(i,1));
    for k=1:numel(prevMinGraphV)
        minV=find(minGraphV==prevMinGraphV(k),1,'first');
        fullV=find(fullGraphV==prevFullGraphV(k),1,'first');
        minGraphVtoFullGraphV(minV)=fullV;
        fullGraphVtoMinGraphV(fullV)=minV;
    end
    
    tmpAllMatch = RecursComputeMatching(minGraphVtoFullGraphV,fullGraphVtoMinGraphV,...
                                        minGraphV,minGraphVSignature,minGraphVDistanceMtr,...
                                        fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr,...
                                        allMatch,1);
    if size(tmpAllMatch,2) >= 1
        break;
    end
end

if size(tmpAllMatch,2)==0
    allMatch={};
else
    allMatch=prevAllMatch(i:size(prevAllMatch,1),:);
end

end



% Recusive function 
function [ allMatch ]= RecursComputeMatching(minGraphVtoFullGraphV,fullGraphVtoMinGraphV,...
                                             minGraphV,minGraphVSignature,minGraphVDistanceMtr,...
                                             fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr,...
                                             allMatch,nbMaxMatching)
global edgeRatioThreshold;
global matchVThreshold;

% Find the next vertex of the minGraph to match
[ nextMinV ] = find(minGraphVtoFullGraphV==0,1,'first');
if isempty(nextMinV)  % If all the vertices have been matched
    % Add this solution to the variable allMatch
    allMatch(size(allMatch,1)+1,1)={fullGraphV(minGraphVtoFullGraphV)};
    return
end


% Compute all possible matching with the newly added vertex and the connected edges
for nextFullV=1:numel(fullGraphV)
    % If the full graph vertex has not been yet matched
    if fullGraphVtoMinGraphV(nextFullV)~=0
        continue;
    end
    % Check if the signature matches minGraphVtoFullGraphV
    if CompVertSimilarity(minGraphVSignature(nextMinV),fullGraphVSignature(nextFullV)) < matchVThreshold
        continue;
    end
    % For all connected vertices to nextV
    goodMatch=true;
    for connectMinV=1:numel(minGraphV)
        % If it has no matching, we skip
        if minGraphVtoFullGraphV(connectMinV)==0
            continue;
        end
        connectFullV=minGraphVtoFullGraphV(connectMinV);
        
        % Check if the edge lengths match
        diffDist=abs(minGraphVDistanceMtr(connectMinV,nextMinV)-...
                     fullGraphVDistanceMtr(connectFullV,nextFullV));
        diffRatio=diffDist*2/(minGraphVDistanceMtr(connectMinV,nextMinV)+...
                              fullGraphVDistanceMtr(connectFullV,nextFullV));
        if diffRatio>edgeRatioThreshold
            goodMatch=false;
            break;            
        end
    end
    if ~goodMatch
        continue
    end
    % We found a matching of nextMinV in the full graph
    nextMinGraphVtoFullGraphV=minGraphVtoFullGraphV;
    nextFullGraphVtoMinGraphV=fullGraphVtoMinGraphV;
    nextMinGraphVtoFullGraphV(nextMinV)=nextFullV;
    nextFullGraphVtoMinGraphV(nextFullV)=nextMinV;
    allMatch=RecursComputeMatching(nextMinGraphVtoFullGraphV,nextFullGraphVtoMinGraphV,...
                                   minGraphV,minGraphVSignature,minGraphVDistanceMtr,...
                                   fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr,...
                                   allMatch,nbMaxMatching);
    % Do not find all the matching, stop when 2 matching are found
    if size(allMatch,2) >= nbMaxMatching
        break;
    end
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add one vertex to the min graph %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ minGraphV,markedV ] = AddVToMinGraph(minGraphV,fullGraphV,fullGraphVDistanceMtr,markedV)

% Extend the minimal graph with one vertex
% The largest geodesic distance to the minimal graph should be the smallest
minDist=realmax;
minDistIdx=0;

for fullV=1:numel(fullGraphV)
    % Find a vertex that does not belong to the minimal graph
    if ~isempty(find(minGraphV==fullGraphV(fullV),1))
        continue;
    end
    
    if markedV(fullV)
        continue;
    end
    % Compute the largest geodesic distance to the vertices of the minGraph
    tmp=0;
    for minV=1:numel(minGraphV)
        % Find the corresponding vertex in the full graph
        fullV2=find(fullGraphV==minGraphV(minV),1);
        if isempty(fullV2)
            continue;
        end
        if fullGraphVDistanceMtr(fullV,fullV2)>tmp
            tmp=fullGraphVDistanceMtr(fullV,fullV2);
        end
    end
    if tmp==0
        continue;
    end
    if tmp<minDist
        minDist=tmp;
        minDistIdx=fullV;
    end;
end

% If no vertex has been found
if minDistIdx==0
    % We return
    return;
end

%disp(['Adding vertex: ',num2str(fullGraphV(minDistIdx))]);

% We add the vertex to the min graph
minGraphV=[minGraphV;fullGraphV(minDistIdx)];
markedV(minDistIdx)=true;

end


% Compute the quality of the matching (difference of edge length and vertex signature)
function [ cost ]= MatchingCost(allMatch,minGraphV,minGraphVSignature,minGraphVDistanceMtr,fullGraphV,fullGraphVSignature,fullGraphVDistanceMtr)

nbMatch=size(allMatch,1);
nbMinGraphV=size(minGraphVSignature,1);

cost=zeros(nbMatch,2);

for i=1:nbMatch
    matchedFullGraphV=cell2mat(allMatch(i,1));
    % Compute the quality of matching for the vertex signature
    vSignatureDiff=0;
    edgeDiff=0;
    for minV=1:nbMinGraphV
        fullV=find(fullGraphV==matchedFullGraphV(minV),1,'first');
        diff=fullGraphVSignature(fullV,:)-minGraphVSignature(minV,:);
        d=diff*diff';
        vSignatureDiff=vSignatureDiff+d;
    end
    % Compute the quality of matching for edge length (geodesic path)
    for minV1=1:nbMinGraphV
        for minV2=minV1+1:nbMinGraphV
            if minGraphVDistanceMtr(minV1,minV2)==0
                continue;
            end
            fullV1=find(fullGraphV==matchedFullGraphV(minV1),1,'first');
            fullV2=find(fullGraphV==matchedFullGraphV(minV2),1,'first');
            diff=minGraphVDistanceMtr(minV1,minV2)-fullGraphVDistanceMtr(fullV1,fullV2);
            d=diff*diff';
            %disp(['edge diff',num2str(d)]);
            edgeDiff=edgeDiff+d;
        end
    end
    cost(i,1)=vSignatureDiff;
    cost(i,2)=edgeDiff;
end
end


function [ minGraphVDistanceMtr ]= CompMinGraphVDistanceMtr(minGraphV,fullGraphV,fullGraphVDistanceMtr)

minGraphVDistanceMtr=zeros(numel(minGraphV),size(fullGraphVDistanceMtr,2));
for minV=1:numel(minGraphV)
    fullV=find(fullGraphV==minGraphV(minV),1,'first');
    minGraphVDistanceMtr(minV,:)=fullGraphVDistanceMtr(fullV,:);
end
end






